Johnson Hsieh
April, 2015 @ DSP A1 coures
Contact me｜ Data Analyst Training Program @ DSP

ubike <- fread("data/ubikeweatherutf8.csv", 
               showProgress = interactive(), data.table = FALSE)
ubike1 <- filter(ubike, sno==1) %>%
  mutate(sbi.range=max.sbi-min.sbi) %>%
  mutate(is.rushhours=cut(hour, breaks=c(0, 8, 10, 17, 20, 24), 
                          labels = c(0,1,0,1,0), right=FALSE)) %>%
  mutate(is.weekday=ifelse(strftime(date, "%u") < 6, 1, 0))

tab1 <- filter(ubike1, is.rushhours==1, is.weekday==1) %>% 
  group_by(tot) %>%
  summarise(min(sbi.range), mean(sbi.range), max(sbi.range))

探索性資料分析 — 以YouBike為例

市府站車輛數基本數據

df1 <- group_by(ubike1, date, hour) %>%
  summarise(rate.sbi=mean(avg.sbi)/tot) %>%
  group_by(hour) %>%
  summarise(rate.sbi=mean(rate.sbi))

thm <- function() {
  theme_gray(base_family = "STHeiti") + 
  theme(text=element_text(size=18))
}

ggplot(df1, aes(x=hour, y=rate.sbi)) +
  geom_bar(stat="identity") + 
  ggtitle("捷運市政府站") +
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank()) + 
  scale_y_continuous(labels = percent)

上圖為捷運市府站每天有車率的變化，大約晚間十點至隔天凌晨兩點間有車率 (當時段平均車輛數 / 總車輛數) 最低，係因此時YouBike公司將車輛回收，於清晨三點左右陸續將車輛補回。由圖可知三點與四點時有車率大幅增加兩次，推測是市府站規模較大需要兩次小貨車補給 (每次約補給30輛車)。在上午七點左右，有車率開始明顯下降，直至上午十點到達低點，即31%。之後有車率略微增加，直至下午一點達到當日次高峰。可以發現，下午三點之後有車率再度明顯下降，直至晚間九點。

場站有車率與晴雨關係

df2 <- filter(ubike, sno==1) %>% 
  mutate(is.rain=rainfall>1) %>%
  mutate(is.rain=factor(is.rain, levels=c(FALSE, TRUE), 
                        labels = c("晴天","雨天"))) %>%
  select(date, hour, tot, avg.sbi, avg.bemp, temp, is.rain) %>%
  group_by(date, hour, is.rain) %>%
  summarise(rate.sbi=mean(avg.sbi)/tot) %>%
  group_by(hour, is.rain) %>%
  summarise(rate.sbi=mean(rate.sbi))

kable(df2[17:22,], format="html", digits=3, align="c")
hour is.rain rate.sbi
8 晴天 0.431
8 雨天 0.658
9 晴天 0.288
9 雨天 0.576
10 晴天 0.254
10 雨天 0.547

長條圖範例

首先用長條圖 (bar chart) 來探索這份報表，當欄位大於二時，將依賴顏色做區隔，一般而言長條圖有以下變型： - Stack plot (堆疊圖) - Dodge plot () - Fill plot (相對堆疊) - Facet panels - Pyramid (金字塔圖)


  • Stack plot (堆疊圖)

    Hint: geom_bar(stat="identity", position="stack")
ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
  geom_bar(stat="identity") + 
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank())


  • Dodge plot

    Hint: geom_bar(stat="identity", position="dodge")
ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
  geom_bar(stat="identity", position="dodge") + 
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank())


  • Fill plot

    Hint: geom_bar(stat="identity", position="fill")
ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
  geom_bar(stat="identity", position="fill") + 
  labs(x="時間", y="相對有車率") +
  thm() +
  theme(legend.title=element_blank())


  • Facet panels in a grid

    Hint: facet_grid(y~.) or facet_grid(.~x)
ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
  geom_bar(stat="identity", position="dodge") + 
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank()) + 
  facet_grid(is.rain~.)

ggplot(df2, aes(x=hour, y=rate.sbi, fill=is.rain)) +
  geom_bar(stat="identity", position="dodge") + 
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank()) + 
  facet_grid(.~is.rain)


  • Pyramid

    Hint: filter(df2, is.rain=="晴天"), and coord_flip()
ggplot(df2, aes(x=hour,y=rate.sbi, fill=is.rain)) + 
  geom_bar(data=filter(df2, is.rain=="晴天"), stat="identity") + 
  geom_bar(aes(y=rate.sbi*(-1)), data=filter(df2, is.rain=="雨天"), 
           stat="identity") + 
  scale_y_continuous(breaks=seq(from=-1, to=1, by=0.1), 
                     labels=abs(seq(-1, 1, 0.1))) + 
  labs(x="時間", y="有車率") +
  theme(legend.title=element_blank()) + 
  coord_flip() + thm()


折線圖範例

由於x軸的單位為時間，所以可以改用折線圖 (line chart) 來做探索，本例就比較晴、雨天有車率的差異而言，line chart 比 bar chart 更具視覺上優勢。


  • Line chart

    由以下折線圖可知，不管在哪個時段，晴天的有車率都低於雨天，而差異最大的時段是午夜十二點前後 (Why?)，最小的時段是凌晨五點前後 (Why?)。
ggplot(df2, aes(x=hour, y=rate.sbi, colour=is.rain, fill=is.rain)) +
  geom_line(size=1) + 
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank())


  • Smoothing (平滑曲線)

    透過統計方法 (loess method)，可以迅速地描繪出平滑的趨勢線 (陰影處為95%信賴區間)。
    Hint: stat_smooth()
ggplot(df2, aes(x=hour, y=rate.sbi, colour=is.rain, fill=is.rain)) +
  #geom_line() + 
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank()) +
  stat_smooth(size=1)


  • Polar coordinates (極座標化)

    另一種呈現方式，又稱雷達圖 (radar chart)，賦予資料面積的概念。
    Hint: coord_polar()
ggplot(df2, aes(x=hour, y=rate.sbi, colour=is.rain, fill=is.rain)) +
  #geom_line() + 
  labs(x="時間", y="有車率") +
  thm() +
  theme(legend.title=element_blank()) +
  stat_smooth(size=1) + coord_polar()

熱點圖範例

熱點圖 (heatmap) 是用顏色深淺呈現數值大小的視覺化。
Hint: geom_tile()

  ggplot(df2, aes(x=hour, y=is.rain, fill=rate.sbi)) + 
  geom_tile() +
  scale_fill_gradient(name="有車率", low="white", high="midnightblue") + 
  labs(x="時間", y="天氣") +
  thm()

平行座標圖

平行座標圖 (Parallel coordinate plot) 多用於呈現多欄位的資料視覺化，強調欄位的順序性，特別適合用在因果關係的陳述。譬如：行業別 -> 是否上DSP課程 -> 職場表現。
Hint: library(GGally) and ggparcoord()

df2 <- mutate(df2, rain=as.numeric(is.rain)-1)
ggparcoord(as.data.frame(df2), columns = c(1,4,3), groupColumn = 2) +
  thm() + theme(legend.title=element_blank())

ggparcoord(data = iris, columns = 1:4, groupColumn = 5,
           title = "Parallel Coordinate Plot for the Iris Data") + thm()

場站 有車率 與 使用率 的探索

df3 <- filter(ubike1, sno==1) %>% 
  mutate(is.rain=rainfall>1) %>%
  mutate(is.rain=factor(is.rain, levels=c(FALSE, TRUE), 
                        labels = c("晴天","雨天"))) %>%
  mutate(is.weekday=strftime(date, "%u")<6) %>%
  mutate(is.weekday=factor(is.weekday, levels=c(FALSE, TRUE),
                           labels=c("平日","假日"))) %>%
  mutate(is.rushhours=cut(hour, breaks=c(0, 4, 7, 24), right=FALSE)) %>%
  group_by(date, is.weekday, is.rushhours, hour, is.rain) %>%
  summarise(rate.sbi=mean(avg.sbi)/tot, rate.used=mean(max.sbi-min.sbi)/tot)

kable(head(df3), format="html", digits=3, align="c")
date is.weekday is.rushhours hour is.rain rate.sbi rate.used
2014-12-08 假日 [7,24) 15 雨天 0.536 0.033
2014-12-08 假日 [7,24) 16 雨天 0.545 0.044
2014-12-08 假日 [7,24) 17 雨天 0.547 0.039
2014-12-08 假日 [7,24) 18 雨天 0.573 0.044
2014-12-08 假日 [7,24) 19 雨天 0.556 0.100
2014-12-08 假日 [7,24) 20 雨天 0.737 0.100

散佈圖

散佈圖 (scatterplot) 是比較兩數值變數最直覺的視覺化，主要是觀察兩變數間是否存在特殊的趨勢、群聚現象、離群值。以下方式常用來輔助散佈圖的探索：

  • 重新定義象限
  • 改變座標軸的尺度 (對數轉換)
  • 加入額外資訊 (大小、顏色、形狀)
  • 各種分群演算法 (K-means)
ggplot(df3, aes(x=rate.used, y=rate.sbi)) + geom_point() + 
  labs(x="使用率", y="有車率") + thm()

  • 定義象限 

ggplot(df3, aes(x=rate.used, y=rate.sbi)) + geom_point() + 
  labs(x="使用率", y="有車率") + thm()+
  geom_vline(xintercept=0.4, lty=2) +
  geom_hline(yintercept=0.4, lty=2)

  • 改變座標尺度 (摩爾定律)


  • 利用顏色分類

    Hint: aes(colour=is.rushhours)
ggplot(df3, aes(x=rate.used, y=rate.sbi)) + 
  labs(x="使用率", y="有車率") + thm() +
  geom_point(aes(colour=is.rushhours), position="jitter")


  • Facet grid, color, and shape

    Hint: aes(colour=is.rushhours, shape=is.weekday)
ggplot(df3, aes(x=rate.used, y=rate.sbi)) + 
  geom_point(aes(colour=is.rushhours, shape=is.weekday), position="jitter") +
  facet_grid(is.weekday~is.rushhours) + 
  labs(x="使用率", y="有車率") +
  thm() +
  theme(legend.title=element_blank())

df4 <- filter(df3, hour>6 & hour < 10)

ggplot(df4, aes(x=rate.used, y=rate.sbi)) + 
  geom_point(aes(colour=paste(as.character(is.weekday),
                              as.character(is.rain), sep="-"))) +
  ggtitle("每日7-9點 YouBike使用狀況") +
  labs(x="使用率", y="有車率") +
  facet_grid(is.rain~is.weekday) + 
  thm() +
  theme(legend.title=element_blank())

與鄰近場站的關係

tmp <- group_by(ubike, sno, sna, sarea, lat, lng) %>% distinct 
dist <- round(distm(x=tmp[, c("lng","lat")])[,1])
df5 <- tmp %>% select(sno, sna, sarea, lat, lng) %>% 
  cbind(dist) %>% arrange(dist) %>% top_n(10, wt = -dist)

研究完單一場站之後，試著探索市府站與鄰近場站的關係，此時需要透過經緯度計算場站與場站之間的距離。透過geosphere套件中的distm函數可以批次計算所有場站之間的兩兩距離，整理得到下表，離捷運市府站最近的場站依序是台北市政府 (438m), 興雅國中 (484m)…。
Hint: library(geosphere), distm, group_by, distinct 

kable(df5, format="html", digits=3, align="c")
sno sna sarea lat lng dist
1 捷運市政府站(3號出口) 信義區 25.041 121.568 0
3 台北市政府 信義區 25.038 121.565 438
5 興雅國中 信義區 25.037 121.569 484
25 永吉松信路口 信義區 25.045 121.572 659
6 世貿二館 信義區 25.035 121.566 718
150 松德公園 信義區 25.037 121.573 734
138 捷運永春站(2號出口) 信義區 25.041 121.575 754
8 世貿三館 信義區 25.035 121.564 759
113 仁愛逸仙路口 信義區 25.038 121.561 763
4 市民廣場 信義區 25.036 121.562 778

地圖應用範例

利用ggmap套件導入google map作為底圖將場站位置標示出來。
Hint: library(ggmap), map <- get_map("Taipei"); ggmap(map), geom_point

df5$is.cityhall <- factor(c(1, rep(0, 9)), levels=1:0)
map <- get_map(location=c(lon=df5$lng[1], lat=df5$lat[1]) , zoom = 15)
ggmap(map) + thm() +
  geom_point(data=df5, aes(x=lng, y=lat, colour=is.cityhall), size=5) + 
  geom_text(data=df5, aes(x=lng, y=lat, label=sna, colour=is.cityhall), 
            position="jitter", vjust=-1, hjust=0.5, size=4, family="STHeiti") + 
  theme(legend.position="none") + scale_color_brewer(palette="Set1")


練習用geom_point(size=tot)來改變場站標示的大小。

df5 <- group_by(tmp) %>% select(sno, tot) %>% 
  right_join(df5, by="sno") %>%
  `[`(c(1, 3, 4, 2, 5, 6, 7, 8)) 
map <- get_map(location=c(lon=df5$lng[1], lat=df5$lat[1]) , 
               maptype = "roadmap", zoom = 15)

ggmap(map) + thm() +
  geom_point(data=df5, aes(x=lng, y=lat, colour=is.cityhall, size=tot)) + 
  theme(legend.position="none") + scale_color_brewer(palette="Set1") +
  scale_size(range = c(3,12))


熱點圖進階應用

觀察鄰近捷運市府站的10個YouBike場站，每一天 有車率 與 使用率的狀況。以有車率為例，透過觀察可以發現{興雅國中, 台北市政府, 市民廣場, 世貿三館, 世貿二館} 時間分佈有相似的狀況，{永吉松信路口, 松德公園, 捷運永春站} 也有相似的情況，而捷運市府站介於兩群之間，仁愛逸仙路口則是一枝獨秀。

  • 使用率的熱點圖可以看到什麼趨勢？
  • 有沒有自動排序的統計方法？
tmp1 <- filter(ubike, sno%in%df5$sno) %>%
  mutate(is.rain=rainfall>1) %>%
  mutate(is.rain=factor(is.rain, levels=c(FALSE, TRUE), 
                        labels = c("晴天","雨天"))) %>%
  mutate(is.weekday=strftime(date, "%u")<6) %>%
  mutate(is.weekday=factor(is.weekday, levels=c(FALSE, TRUE),
                           labels=c("平日","假日"))) %>%
  mutate(is.rushhours=cut(hour, breaks=c(0, 4, 7, 24), right=FALSE)) %>%
  group_by(date, sno, sna, is.weekday, is.rushhours, is.rain, hour, tot) %>%
  summarise(rate.sbi=mean(avg.sbi)/tot, rate.used=mean(max.sbi-min.sbi)/tot)

df6 <- tmp1 %>%
  filter(is.weekday=="平日", is.rain=="晴天") %>%
  group_by(sno, sna, sna, hour) %>%
  summarise(rate.sbi=mean(rate.sbi), rate.used=mean(rate.used))

ggplot(df6, aes(x=hour, y=sna, fill=rate.sbi)) + geom_tile() + thm() + 
  theme(legend.position="bottom") + 
  scale_fill_gradient(name="有車率", low="white", high="lawngreen") + 
  labs(x="時間", y="") +
  theme(axis.text = element_text(size = 13, color="darkgreen"))

ggplot(df6, aes(x=hour, y=sna, fill=rate.used)) + geom_tile() + thm() + 
  theme(legend.position="bottom") + 
  scale_fill_gradient(name="使用率", low="white", high="Navy") + 
  labs(x="時間", y="") +
  theme(axis.text = element_text(size = 13, color="darkblue"))


  • heatmap 排序

    當heatmap的x軸或y軸為類別變數時，可以經由階層分群法 (hierarchical clustering) 做行或列的排序。

  • 首先我們需要一個 場站對時間 (sna ~ hour) 的有車率 (rate.sbi) 矩陣 (dcast)
  • 使用階層分群演算法 (hclust)
  • 畫出分群樹狀圖 (ggdendrogram)

  • 取得排序 (order)

Hint: library(reshape2), library(ggdendro)

dat <- dcast(df6, sna~hour, value.var="rate.sbi")
rownames(dat) <- dat[,1]
dat <- dat[,-1]
kable(dat[,8:13], format = "html", digits = 3, row.names = TRUE)
7 8 9 10 11 12
捷運市政府站(3號出口) 0.589 0.553 0.504 0.463 0.479 0.484
捷運永春站(2號出口) 0.276 0.229 0.188 0.192 0.160 0.156
仁愛逸仙路口 0.228 0.184 0.157 0.151 0.137 0.192
世貿二館 0.148 0.216 0.372 0.522 0.760 0.767
世貿三館 0.146 0.197 0.270 0.397 0.469 0.465
市民廣場 0.204 0.257 0.316 0.389 0.487 0.510
松德公園 0.267 0.218 0.156 0.095 0.059 0.039
台北市政府 0.251 0.242 0.297 0.386 0.497 0.625
興雅國中 0.192 0.287 0.418 0.543 0.627 0.686
永吉松信路口 0.280 0.245 0.166 0.092 0.058 0.049 
hc.sna <- hclust(dist(dat))
ggdendrogram(hc.sna, rotate = TRUE) + thm() + labs(x="", y="")

# hc.sna$order
sna.order <- data.frame(order=1:10, sna=hc.sna$labels[hc.sna$order])
kable(sna.order, format = "html")
order sna
1 仁愛逸仙路口
2 捷運永春站(2號出口)
3 松德公園
4 永吉松信路口
5 捷運市政府站(3號出口)
6 世貿三館
7 市民廣場
8 世貿二館
9 台北市政府
10 興雅國中 
df7 <- df6
df7$sna <- factor(df7$sna, levels=(sna.order[,2]))
ggplot(df7, aes(x=hour, y=sna, fill=rate.sbi)) + geom_tile() + thm() + 
  theme(legend.position="bottom") + 
  scale_fill_gradient(name="有車率", low="white", high="lawngreen") + 
  labs(x="時間", y="") +
  theme(axis.text = element_text(size = 13, color="darkgreen"))


  • 對時間做排序 

hc.hour <- hclust(dist(t(dat)))
ggdendrogram(hc.hour) + thm() + labs(x="", y="")

hour.order <- data.frame(order=1:24, sna=hc.hour$labels[hc.hour$order])
# kable(hour.order, format = "html")

df7$hour <- factor(df7$hour, levels=(hour.order[,2]))

ggplot(df7, aes(x=hour, y=sna, fill=rate.sbi)) + geom_tile() + thm() + 
  theme(legend.position="bottom") + 
  scale_fill_gradient(name="有車率", low="white", high="lawngreen") + 
  labs(x="時間", y="") +
  theme(axis.text = element_text(size = 13, color="darkgreen"))


  • 試著對 使用率 進行排序 

dat <- dcast(df6, sna~hour, value.var="rate.used")
rownames(dat) <- dat[,1]
dat <- dat[,-1]
hc.sna <- hclust(dist(dat))
hc.hour <- hclust(dist(t(dat)))
df8 <- df6
df8$sna <- factor(df8$sna, levels = hc.sna$labels[hc.sna$order])
df8$hour <- factor(df8$hour, levels = hc.hour$labels[hc.hour$order])
ggplot(df8, aes(x=hour, y=sna, fill=rate.used)) + geom_tile() + thm() + 
  theme(legend.position="bottom") + 
  scale_fill_gradient(name="使用率", low="white", high="Navy") + 
  labs(x="時間", y="") +
  theme(axis.text = element_text(size = 13, color="darkblue"))


平行座標圖進階應用

平行座標圖常用來展示不同群組在諸多變數間的差異性，當群組分類方式未知時，可以利用機器學習 (machine learning) 中的非監督式學習 (unsupervised learning)，幫資料做分群。分群之後再藉由平行座標圖來呈現資料的脈絡。

  • 選擇 平日, 晴天, 7-21時鄰近市府站的資料進行分析
  • 以場站大小 (tot)、有車率 (rate.sbi)、使用率 (rate.used) 三個變數做分群
  • 使用K-means演算法分3群
  • 將分群結果視作新的變數畫平行座標圖

tmp2 <- filter(tmp1, is.weekday=="平日", is.rain=="晴天", hour>6 & hour<22) %>%
  group_by(sno, sna, tot) %>%
  summarise(rate.sbi=mean(rate.sbi), rate.used=mean(rate.used))
km <- kmeans(tmp2[,3:5], 3)
km
K-means clustering with 3 clusters of sizes 1, 4, 5

Cluster means:
    tot  rate.sbi rate.used
1 180.0 0.3260354 0.1156761
2  65.0 0.3684294 0.1783035
3  35.2 0.2290803 0.2471085

Clustering vector:
 [1] 1 3 2 2 2 2 3 3 3 3

Within cluster sum of squares by cluster:
[1]   0.0000 300.0137  92.8441
 (between_SS / total_SS =  97.8 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"    
[5] "tot.withinss" "betweenss"    "size"         "iter"        
[9] "ifault"      
df9 <- group_by(tmp2) %>% 
  transmute(sna, tot, rate.sbi, rate.used,
            group=factor(km$cluster))  %>% 
  arrange(group)

ggparcoord(as.data.frame(df9), columns = c(1,2,3,4), groupColumn = 5,
           scale="uniminmax") + 
  geom_line(size=1) + thm() + theme(legend.title=element_blank()) +
  scale_x_discrete(labels=c("場站","總停車格","有車率","使用率")) +
  labs(x="", y="")

kable(df9, format="html", digits=3, align="c")
sna tot rate.sbi rate.used group
捷運市政府站(3號出口) 180 0.326 0.116 1
市民廣場 60 0.367 0.137 2
興雅國中 60 0.384 0.163 2
世貿二館 80 0.410 0.257 2
世貿三館 60 0.313 0.156 2
台北市政府 40 0.361 0.243 3
永吉松信路口 30 0.200 0.266 3
仁愛逸仙路口 38 0.180 0.192 3
捷運永春站(2號出口) 30 0.234 0.356 3
松德公園 38 0.170 0.179 3